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An important aspect of molecular fluids is the relation between orientation and translation parts 
of the two-particle correlations. Especially the detailed knowledge of the influence of orientation 
correlations is needed to explain and calculate in detail the occurrence of a nematic phase. 
The simplest model system which shows both orientation and translation correlations is a system 
of hard ellipsoids. We investigate an isotropic fluid formed of hard ellipsoids with Percus- Yevick 
theory. 

Solving the Percus- Yevick equations self-consistently in the high density regime gives a clear crite- 
rion for a nematic instability. We calculate in detail the equilibrium phase diagram for a fluid of 
hard ellipsoids of revolution. Our results compare well with Monte Carlo Simulations and density 
functional theory. 
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I. INTRODUCTION 

Isotropic simple liquids formed of atomic systems with 
rotational symmetry are well understood. If the two par- 
ticle correlation is given by a hard sphere interaction an 
integral equation like the Percus- Yevick (PY) 0] closure 
relation can be solved analytically. For the liquid phase 
the PY equation gives good results even in the dense liq- 
uid regime (up to a packing fraction (j) < 0.49) above 
which the equilibrium state is crystalline which PY fails 
to obtain. The packing fraction (p is defined as the rela- 
tion between the number density p and the volume of the 
particles 4> = ^p<J^ with a usually set to one. For hard 
spheres w 0.64 corresponds to random closed packing 
and (j) = V2tt/6 ~ 0.74 to packing in a fee lattice. 

Molecular systems have usually complicated potentials 
which are modeled e.g. by Lennard-Jones potentials of 
each atom in a molecule. A very basic feature of molec- 
ular systems is the existence of orientation degrees of 
freedom which interplay in non-trivial ways with trans- 
lational degrees of freedom. The simplest model system 
which allows to study this interplay is a system of rota- 
tional symmetric hard ellipsoids. The equilibrium phase 
diagram is now - compared to hard spheres - enriched 
by an additional variable, the aspect ratio Xq of the el- 
lipsoids, which is defined as the ratio between the major 
axis, a, and the minor axis, b, Xq = f • As a function of 
the aspect ratio or density a fluid of hard ellipsoids can 
now also undergo an isotropic to nematic (I-N) transi- 
tion. This is expected from the Onsager solution of hard 
sphcrocylinders |^ , has been found by computer simula- 
tions and also by density functional theory 

PY theory deals with the isotropic phase. It is not able 
to describe a phase transition. However it is known from 
the hyper-netted-chain (HNC) closure relation for hard 



ellipsoids |5| and for dipolar hard spheres ||6[ that it is 
in principle possible to identify a precursor phenomenon 
of a phase transition in an integral equation. Such a 
precursor phenomenon is not known for the PY closure 
relation (see also 0). 

In this paper however we show that the PY closure 
relation leads not only to a clear indication of a nematic 
instability but enables us also to calculate the equilibrium 
phase diagram of hard ellipsoids. 

II. INTEGRAL EQUATIONS 

A fundamental relation which all closure relations for 
integral equations are based upon is the Ornstein-Zernike 
(OZ) equation @ 

/i(ri,rJi,r2,fi2) = c(ri,fii,r2,fi2) 

+ p(c{ri,ni,r3,n3)h{r3,n3,r2,n2)) . (1) 

It provides the relation between the total correlation 
function h{ri, r2, £^2) and the direct correlation func- 
tion c(ri, ill, r2, fi2)- The product is given by 

(■■■)„,,, = IS?/""'/""- ■ 

r; is the position of the center of mass of the ellip- 
soid i and fli is the orientation of this ellipsoid repre- 
sented by the Euler angles <I>i , 9i (the third Euler angle 
X is not needed due to the symmetry of the ellipsoids). 
Due to translational invariance the functions depend on 
Tij = Ti — Vj only and due to the isotropy (in the liquid) 
the correlation functions only depend on r = |r.y | once 
a specific coordinate system has been chosen. In the fol- 
lowing subsections we transform the OZ equations and 
the PY closure relation (eq. (|lj)) in such a way that we 
can use it for a numerical treatment. With most of the 
definitions we follow the book of Gray and Gubbins [0] . 
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A. Spherical harmonics 

An obvious orthogonal basis set to expand the angu- 
lar dependence of the correlation functions are spher- 
ical harmonics. The transformed correlation function 
F e {c, h, ...} is given by: 

F{li,l2,m;r) 

= [ dn^ /d^!2i^(^^l,^^2,r)y^7(^!l)y^^*(^!2) (3) 



In the final step one goes from the rotational invariant 
representation to the q-frame and one gets a representa- 
tion of F. 

F{li,l2,m;q) = 



E \I^^^TZ^ ^('i' ^2, 1; q) C{k,h, I; m, -m, 0) (7) 



47r 



Note that our transformation differs by a factor of (—1)™ 
and by a factor of id"'^) fj-om the definition used in the 
book of Gray et al. The latter of these two factors 
gives us (in q-frame) only real elements of the correlation 
functions. For eq. (j^) the r-frame was used. This means 
the z-axis of the coordinate system was chosen along the 
axis connecting the two particles which are correlated [|| . 
Therefore we only need to deal with one index m. In q- 
space we use, after Fourier transformation the laboratory 
fixed frame, the g-frame |^,^ where now all g-dependent 
correlation functions are diagonal in m. Within the com- 
plete set of spherical harmonics the OZ-equation can be 
rewritten: 

h{li,l2, m; r) = c(Zi, ^2, "i; r) 

+ J-'^ dri c{li,l,m;ri)h{l,l2,m;r - ri) (4) 

This equation relates the total correlation func- 
tion h(li,l2,m;r) with the direct correlation function 
c(li, I2, rn] r). The total correlation function has two con- 
tributions, a direct one which results from direct correla- 
tions and is just c(Zi,Z2,w;r) plus an indirect contribu- 
tion which averages over possible interactions mediated 
by another particle in an indirect way. 



B. Transformation into q-space 

Due to the expansion in spherical harmonics a Fourier 
transform cannot be performed as usual. First one has to 
find a representation of F which is invariant with respect 
to rotation. 



F{hj2,l;r) 



47r 



(2^ + 1) 



F{li,l2,m;r) C{li,l2.,l]m,-m,{)) (5) 



where C(/i, I2, 1] mi, 7712, m) are the the Clebsch-Gordon 
coefficients. The next step is the Hankel transformation 
which uses the Rayleigh expansion to transform F from 
r-space to q-space. This involves spherical Bessel func- 
tions ji (qr) due to the expansion of e*^*" within the basis 
of spherical harmonics. 

/"OO 

F{h,l2,l;q) = M-^y jiiqr) F{h,l2,l;r) (6) 
JO 



Therefore the equations (||) - (0) transform a two particle 
correlation function given in real space and r-frame into 
a function in g-space and g-frame. 



C. The Ornstein Zernike equations in q-space 

Applying eqs. (||) - (|7|) to the OZ equation one can 
rewrite eq. (||): 



Hh,l2,m;q) = c(/i, ^2, "i; g) 
_P_ 
47r 



+ j- E ^7 ™; q)Hl, h, m; q) 



(8) 



This can be written as a matrix equation for each m and 
q value: 

k{m; q) = g{m; q) + ■£^g{m; q)k{m; q) (9) 

c,h are symmetric matrices with indices li,l2- 

For the input into our numerical calculation we define 
an auxiliary correlation function y in the usual way [^: 

y{m;q)=h.{m;q)-Q{m;q) (10) 
Using this auxiliary function the OZ equation rewrites: 

(1 - <?)) gim; 9) = ^ 9)] ^ (11) 

This is a linear system of equations which determines 
y{li, I2, m; q) if c(/i, ^2j m] q) is known. 

III. PERCUS-YEVICK CLOSURE RELATION 



In a formal way one can define a product between two 
correlation functions c — a * h |L0| |. In r-frame this 
product reads: 



z{hj2,m;r) = ^ E 



'(2/; + l)(2^^ + l)(2/;' + l)(2/^' + l) 
(2Zi + l)(2;2 + l) 



C{l[,l'l,h;0,0,0) C{l'2,lll2;0,0,0) 

C(l[,l'i,h\m',m",m) C{l2,l2,l2; ~m' , -m" , -m) 



in' ,rn" 

a{l[,l'2,m';r) b{l'l , I2 , m" ; r) 



(12) 
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The Clebsch-Grodan coefScients enter into the equation 
due to spatial rotations which have to be performed. 
The PY cfosure relation can now be expressed as 



b * g 



(13) 



were g is the pair correlation function and via 5=1 — 
e^" the pair potential u enters into the equation. For 
the purpose of solving the PY equation numerically it is 
better to rewrite eq. ( |l3| ) as a function of the auxiliary 
function y. 



f * y 

were / is the Mayer function 



(14) 



f{VLi,VL2,r) = e'^ "(f2i,f22,r) _ ^ ^^^^ 

The matrix elements of / in the basis set of spherical 
harmonics have to be computed using eq. (^). 

This equation ( [1^ ) determines the direct correlation func- 
tion c if the auxiliary function y and the Mayer function 
/ are known. 



(i) Due to the head-tail symmetry of the ellip- 
soids all matrix elements of a correlation function 
F{li, I2, m,u) u £ {r, q} with li odd are zero. 

(ii) All elements of F are real both in r-space and q- 
space. Using the definition of eq. (||) this is even 
valid for all linear molecules. 

(iii) Therefore there is an additional symmetry 
F{li,l2,m,u) = F{li,l2,~m,u) 

(iv) Also the / occurring in the rotational invariants (eq. 
(H)) can only have even values following from li + 
I2+I =even which results from inversion symmetry. 

For a proof of (ii) and (iii) the reader may consult ||l^ and 
for (i) and (iv) one might look up |^ . There is one further 
feature we want to point out: For small argument (r or 
k), all non diagonal elements (h h) have to vanish, at 
least in the isotropic phase: 



lim F{li, I2, m,u) = , if 7^ I2 



(19) 



This follows from the transformation of F under rotations 
R: 



A. The pair potential 

In order to determine the matrix elements of the Mayer 
function 



lim F(li, I2, m, Ru) 

= ^™o 51 Dt:UR)D'ri^,,„r{R)F{hM.m^.rn2.u) (20) 



f{n^,n2,r) = 



for D{ni,n2,r) < r 
-1 for D{ni,n2,r)>r 



(16) 



we use the well known approximation of Berne and 
Pechukas |ll], were D depends on the relative orienta- 
tion of the two ellipsoids: 



Dini,n2,r) 



1 / (cos 6*1 cos 6*2)^ 



1 + X (6162) 
(cos 9i — cos 6*2)^ 
1 - X (eie2) 



-1/2 



(17) 



are unit vectors along the symmetry axis of an ellipsoid 
on position i. This approximation models the interaction 
between two ellipsoids by the overlap of Gaussians. The 
value of X is related to the aspect ratio of the ellipsoids 
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(18) 



Note that Z?(rii, 172, ?") is not invariant under Xo ^ 
which implies x ~^ ~X- 



This equation has to be valid for all R which results in 
Sii,i2^rni.m2- This Can be seen by integrating both sites 
of the above equation and by making use of the unitarity 
if the rotation matrices: 

/dR lim F{li, I2, m, Ru) 
= 8TT^F{h,l2,m,0) 

= j dRDl\*„XR)D'A2,miR)Fih,l2,m,,m2,0) 



E 



8n' 



2h + l 



Sh,l2Smi,m2F{ll,l2, mi, TO2, 0) (21) 



Therefore F has to vanish for small u for all non-diagonal 
elements {h ^ h)- Further the value of the diagonal 
elements of F{li, li,m,u ^ 0) with {2li + 1) different m 
can (for short ranged potentials) not depend on m. This 
symmetry can clearly be seen from fig. |^c and from fig. 



C. Calculation procedure 



B. Symmetries of the solution 

Due to the symmetries of the ellipsoid there are certain 
simplifications in the calculation which can be applied. 



In order to obtain a numerical solution of the equations 
above the following steps have to be performed: 
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1. ) The matrix elements /(Zi,/2,w,r) of the Mayer 

function have to be computed using eq. (||) and 
eqs. (16) - (HI). For our calculation we used 100 
points in the range Min(l,Xo) < r < Max(l,Xo) 
where r is given in units of the major axis a of the 
ellipsoids. 

2. ) An initial guess for c{li^l2,m,r) has to be made 

and a grid for r of N^ax — 400 points in the range 
< r < 10 was chosen. 

3. ) The iteration begins by using eqs. (^ - (0) to 

go from c{li,l2,m,r) in r-space and r-frame to 
€{11,12,™, q) in q-space and q-frame using the 
Rayleigh transformation for the direct correlation 
function c{li,l2,'m,r). It turned out to be crucial 
to use an analytic expansion of the spherical Bessel 
functions ji{x) for small argument x. 



4.) The auxiliary function y{li,l2,m,q) has to be cal- 
culated using the OZ eqs. in the form of eq. (11). 
As a q-space grid we used 400 points in the range 
< q < 50 were q is measured in units of [^] . 



5. ) Using eqs. (0) - (||) we get the function 

y(li,l2, m, r) in r-space and r-frame. 

6. ) With the help of the PY eq. (|l2|) one can obtain 

the next iteration for c(/i, I2, m, r). 

7. ) The steps 3.) to 6.) have to be iterated until a 

fix point of the equations has been reached with a 
given accuracy. In this way a self-consistent solu- 
tion can be found. 

As a test for self-consistency we choose the mean square 
deviation of c between two steps of iteration 



IV. RESULTS FROM PY 

The virial expansion of hard ellipsoids of revolution is 
symmetric with respect to Xq — > 1 / Xq up to the second 
order in density. This approximate symmetry is violated 
for higher densities. 
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/ J2 {c^P+^Kluh,m,r^)-dr>){h,l2,m,r,,)f (22) 

where the summation indices were in the regions h,l2 G 
{0, 2, Imax}, m e {0, 1, 2, 3, ...jmrnax} and at the end 
e was typically chosen to be smaller then 2 • 10~^ as a 
condition for convergence. 

In this way one can obtain a stable self-consistent solu- 
tion for the correlation functions. This has already been 
done for a fluid of ellipsoids in refs. ||5|,|l3|-Jl5| and also in 
ref. ||l^ for a single ellipsoid in a fluid of hard spheres. 
In this work we have extended the calculation to a much 
higher density regime than it has been done in previous 
works. In order to reach the high densities we were forced 
to restrict the maximum numbers for Nmax was set to 400 



and I 



and rrimax to 4. The value of / in the rotational 



invariants was restricted to I G {0, 2, 4, ...,2lmax}- 
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FIG. 1. In order to demonstrate the breakdown of the 

approximate symmetry between prolate and oblate ellipsoids 

(Xo — ► 1/^0) we plotted for Xq = 2.0 and Xq = 0.5 

elements of S{li,l2,'m,q), (a) h — I2 = m = 0, (b) 

h — I2 — 2,m = 0, (c) h — 2,l2 = m = The q-axis 

1/3 

was scaled by a factor of Xq . 



4 



This is shown in fig. |l| where we have plotted three 
matrix elements 3(0,0,0, q), S{2, 2,0, q) and S{2,0,0,q) 
of the static structure factor S{li,l2,m,q) for — 
0.62 which is related to the total correlation function 
h{li,l2,m,q) by 

S{li,l2,m,q) = Si-^j^ + -^h{li,l2,m,q) 

4:71 

Sim, q) = 1 + -r-Zi(™, q) 
— — At:— 

^Kg)- (i-£g(m,g))"' (23) 

The q-axis has been stretched by a factor of It 
can be clearly seen that a symmetry Xq — > 1 / Xq is not 
exactly valid at such high densities. 

A. Nematic instability 

Close to the nematic instability the matrix element 
S{2,2,m,q) of the static structure factor develops a di- 
vergence at q — > 0. This was already discussed in for 
results based on the HNC (hyper-netted chain) closure 
relation. 
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FIG. 2. The static structure factors of two systems of el- 
lipsoids close to and far away from the nematic instability 
are shown by comparing the S{l\,l2,m,q) components for 
(j) = 0.55 for Xo = 1.3 (far away from the nematic insta- 
bility) and for Xo = 2.5 (close to the nematic instability). In 
part (a) the S{Q, 0, 0, g) is plotted in part (b) 5(2, 2, 0, q) and 
in part (c) 5(2, 0, 0, q). Note that the 5(2, 0, 0, q) components 
vanish at g = due to the symmetries. 

In fig. such a precursor of a divergence is seen where 
for (j) = 0.55 two systems — 1.3 (far from the nematic 
phase) and Xq = 2.5 (close to the nematic instability) are 
compared. For Xq = 1.3 the S{0,0,0,q) matrix element 
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dominates as can be seen from fig. Note also from fig. 
§ that the S'(0, 0,0,9) components and the S{2,2,0,q) 
seem to change their role when going from Xq — 1.3 to 
Xq = 2.5. The static structure for Xq — 1.3 is dominated 
by the S{0, 0, 0, q) component while the oricntational cor- 
relator S'(2, 2, 0, q) is small. This behaviour is reversed 
when looking at Xq — 2.5. Here the static structure is 
dominated by the oricntational correlations S{2,2,0,q) 
while the center of mass S{0, 0, 0, q) only shows a weak 
structure. Due to the fact that all the non-diagonal ele- 
ments of 5 at g ^ vanish (as can be seen in fig. |^c for 
S{2, 0, 0, q)) only the c(2, 2, to, 0) component governs the 
nematic instability and we get as a condition for such an 
instability 



hm ( 1 



47r 



c(2, 2, m, q) 







(24) 



This expression is the inverse of the Kerr constant K for 
non dipolar potentials. 



5-0 



P 



X-^ = lim(l-^c(2,2,TO,g) 



(25) 




FIG. 3. In (a) c(2, 2, m, q) for q — > and for different den- 
sities is plotted in (a). The curves for <j) = 0.53, 0.54, 0.55 were 
obtained by a solution of the PY equations whereas the curves 
for (j) = 0.57, 0.59 were obtained by applying a quadratic ex- 
trapolation to higher densities to c(2,2,0, g). The important 
part for the nematic instability is plotted in (b) where the 
function 1 — ^c(2,2,0, g) is drawn which becomes at q — 
the inverse of 5(2, 2, 0, g = 0). 

A detailed analysis of the q ^ behavior therefore 
gives us a condition for the nematic instability as it is also 
discussed in a similar way in |^] where dipolar fluids in the 
HNC approximation have been considered. The instabil- 
ity is demonstrated in fig. ||. In part (a) for Xq = 2.5 
the function c(2,2,0, 5) is plotted for different densities 
close to the nematic instability. The first three densities 
(j) = 0.53, 0.54, 0.55 were the highest ones we could reach 
with the numerical solution of the PY equation and the 
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two higher densities (p = 0.57, 0.59 are quadratic extrap- 
olations. In fig. ^ this is shown in greater detail where 
1 — ^c(2, 2, 0, q) is plotted which it becomes the inverse 
of S'(2,2,0,0) at g = 0. For Xq = 2.5 The critical den- 
sity where S{2, 2, 0, 0) diverges is (p — 0.593, according to 
such a quadratic extrapolation. 



B. Equilibrium phase diagram 




FIG. 4. Isotropic to nematic instability as it arises from a 
Percus-Yevick calculation. The highest densities considered 
are plotted with triangles. From there an extrapolation was 
done to determine the instability by c(2, 2, 0, g ^ 0) = An/ p. 
This was done with a linear (circles) and a quadratic (squares) 
interpolation. For comparison we plotted the I-N-transition 
which arises from density functional theory Also the re- 
sults from ref. |^ are plotted with black crosses. 

Using the above condition (eq. (p4|)) we get the phase 
diagram for the hard ellipsoids as it arises from the PY 
approximation as shown in fig. ^ For all densities con- 
sidered there results a clear indication for a nematic in- 
stability. This phase boundary still depends on the way 
the extrapolation to higher densities is done (here linear 
or quadratic) but is in good agreement with other works. 
For example density functional theory of Groh and Di- 
etrich [Q are in a reasonable good agreement with our 
results. However due to their approximation the Xq to 
1/Xo symmetry is exact which is however clearly bro- 
ken in the PY result. Also the results of Frenkel et al. 
1^ are along the same line. It also seems to be that 
density-functional theory shows a better agreement with 
PY theory for prolate ellipsoids that for oblate ones. 



V. CONCLUSION 

Orientational degrees of freedom in molecular systems 
can drive a phase transition into an orientationally or- 
dered nematic liquid crystal phase. In principle integral 
equations have the ability to describe a precursor phe- 
nomenon of such an orientational transition. Until now 
this is well known for e.g. hyper-netted chain (HNC) the- 
ory. In this work we demonstrate that also Percus-Yevick 
theory gives a clear precursor of the nematic phase. We 
therefore were able to calculate the equilibrium phase di- 
agram of hard ellipsoids of revolution. 
The obtained phase diagram is in good agreement with 
density functional theory jj] and Monte Carlo simula- 
tions §]. 
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